# Estimate models that examine how margin of victory, opinion poll prediction, and surprising results
# condition the impact of election outcomes on voter emotions. See Appendix A.




#############
# load data #
#############

library(brms)
library(ggplot2)
library(ggpubr)
library(stargazer)

load("PSRM Replication Files/PresidentialElectionsPosts.RData")

colnames(president_df)

all(c("poll_predict_std", "marginvictory_std", "surprise") %in% colnames(president_df)) # TRUE

length(unique(president_df$country)) # 26
length(unique(president_df$electionname)) # 29
length(unique(president_df$party_country)) # 52

# focus on 15 days before and after the election
day_range <- 15

sub <- president_df[which(president_df$daysinceelection >= -1*day_range
                          & president_df$daysinceelection <= day_range),]

# winners and losers
winners <- sub[which(sub$win == 1),]
losers <- sub[which(sub$win != 1),]

length(unique(winners$electionname[which(!is.na(winners$polarization))])) # 27
length(unique(winners$electionname[which(!is.na(winners$economic_polarization))])) # 27
length(unique(winners$electionname[which(!is.na(winners$social_polarization))])) # 27
length(unique(winners$electionname[which(!is.na(winners$populist_present))])) # 27
length(unique(winners$electionname[which(!is.na(winners$populist_dummy))])) # 27

length(unique(losers$electionname[which(!is.na(losers$polarization))])) # 27
length(unique(losers$electionname[which(!is.na(losers$economic_polarization))])) # 27
length(unique(losers$electionname[which(!is.na(losers$social_polarization))])) # 27
length(unique(losers$electionname[which(!is.na(losers$populist_present))])) # 27
length(unique(losers$electionname[which(!is.na(losers$populist_dummy))])) # 27




##################
# winners + love #
##################

winner_love_poll <- brm(loveprop ~ post*poll_predict_std 
                        #+ post*populist_present + post*polarization
                        + populist_present + polarization
                        + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                        + Xl + Xr 
                        + (1 + Xl + Xr | party_election),
                        data=winners,
                        warmup=1000, iter=2000, chains=3, seed=1234)
#plot(winner_love_poll)
#summary(winner_love_poll)

winner_love_margin <- brm(loveprop ~ post*marginvictory_std
                          #+ post*populist_present + post*polarization
                          + populist_present + polarization
                          + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                          + Xl + Xr 
                          + (1 + Xl + Xr | party_election),
                          data=winners,
                          warmup=1000, iter=2000, chains=3, seed=1234)
#plot(winner_love_margin)
#summary(winner_love_margin)

winner_love_surprise <- brm(loveprop ~ post*surprise 
                            #+ post*populist_present + post*polarization
                            + populist_present + polarization
                            + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                            + Xl + Xr 
                            + (1 + Xl + Xr | party_election),
                            data=winners,
                            warmup=1000, iter=2000, chains=3, seed=1234)
#plot(winner_love_surprise)
#summary(winner_love_surprise)




###################
# winners + angry #
###################

winner_angry_poll <- brm(angryprop ~ post*poll_predict_std 
                         #+ post*populist_present + post*polarization
                         + populist_present + polarization
                         + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                         + Xl + Xr 
                         + (1 + Xl + Xr | party_election),
                         data=winners,
                         warmup=1000, iter=2000, chains=3, seed=1234)
#plot(winner_angry_poll)
#summary(winner_angry_poll)

winner_angry_margin <- brm(angryprop ~ post*marginvictory_std 
                           #+ post*populist_present + post*polarization
                           + populist_present + polarization
                           + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                           + Xl + Xr 
                           + (1 + Xl + Xr | party_election),
                           data=winners,
                           warmup=1000, iter=2000, chains=3, seed=1234)
#plot(winner_angry_margin)
#summary(winner_angry_margin)

winner_angry_surprise <- brm(angryprop ~ post*surprise
                             #+ post*populist_present + post*polarization
                             + populist_present + polarization
                             + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                             + Xl + Xr 
                             + (1 + Xl + Xr | party_election),
                             data=winners,
                             warmup=1000, iter=2000, chains=3, seed=1234)
#plot(winner_angry_surprise)
#summary(winner_angry_surprise)




#################
# losers + love #
#################

loser_love_poll <- brm(loveprop ~ post*poll_predict_std  
                       #+ post*populist_present + post*polarization
                       + populist_present + polarization
                       + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                       + Xl + Xr 
                       + (1 + Xl + Xr | party_election),
                       data=losers,
                       warmup=1000, iter=2000, chains=3, seed=1234)
#plot(loser_love_poll)
#summary(loser_love_poll)

loser_love_margin <- brm(loveprop ~ post*marginvictory_std  
                         #+ post*populist_present + post*polarization
                         + populist_present + polarization
                         + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                         + Xl + Xr 
                         + (1 + Xl + Xr | party_election),
                         data=losers,
                         warmup=1000, iter=2000, chains=3, seed=1234)
#plot(loser_love_margin)
#summary(loser_love_margin)

loser_love_surprise <- brm(loveprop ~ post*surprise 
                           #+ post*populist_present + post*polarization
                           + populist_present + polarization
                           + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                           + Xl + Xr 
                           + (1 + Xl + Xr | party_election),
                           data=losers,
                           warmup=1000, iter=2000, chains=3, seed=1234)
#plot(loser_love_surprise)
#summary(loser_love_surprise)




##################
# losers + angry #
##################

loser_angry_poll <- brm(angryprop ~ post*poll_predict_std
                        #+ post*populist_present + post*polarization
                        + populist_present + polarization
                        + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                        + Xl + Xr 
                        + (1 + Xl + Xr | party_election),
                        data=losers,
                        warmup=1000, iter=2000, chains=3, seed=1234)
#plot(loser_angry_poll)
#summary(loser_angry_poll)

loser_angry_margin <- brm(angryprop ~ post*marginvictory_std
                          #+ post*populist_present + post*polarization
                          + populist_present + polarization
                          + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                          + Xl + Xr 
                          + (1 + Xl + Xr | party_election),
                          data=losers,
                          warmup=1000, iter=2000, chains=3, seed=1234)
#plot(loser_angry_margin)
#summary(loser_angry_margin)

loser_angry_surprise <- brm(angryprop ~ post*surprise 
                            #+ post*populist_present + post*polarization
                            + populist_present + polarization
                            + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                            + Xl + Xr 
                            + (1 + Xl + Xr | party_election),
                            data=losers,
                            warmup=1000, iter=2000, chains=3, seed=1234)
#plot(loser_angry_surprise)
#summary(loser_angry_surprise)


#save(winner_love_surprise, winner_love_poll, winner_love_margin,
#     winner_angry_surprise, winner_angry_poll, winner_angry_margin,
#     loser_love_surprise, loser_love_poll, loser_love_margin,
#     loser_angry_surprise, loser_angry_poll, loser_angry_margin,
#     file="Stan_Surprise.RData")




##################
# result summary #
##################

load("PSRM Replication Files/Stan_Surprise.RData")

summary(winner_love_margin)
summary(winner_angry_margin)
summary(loser_love_margin)
summary(loser_angry_margin)

summary(winner_love_poll)
summary(winner_angry_poll)
summary(loser_love_poll)
summary(loser_angry_poll)

summary(winner_love_surprise)
summary(winner_angry_surprise)
summary(loser_love_surprise)
summary(loser_angry_surprise)


# plot interaction terms
getStdPostCI <- function(model, var1, var2, intr, label){
  coeftab <- summary(model)[[14]][,c(1,3,4)]
  coeftab <- coeftab[intr,]
  
  temp <- data.frame(label=label,
                     est=coeftab[1],
                     lwr=coeftab[2],
                     upr=coeftab[3])
  colnames(temp) <- c("label", "est", "lwr", "upr")
  return(temp)
}

citab <- rbind(getStdPostCI(winner_love_margin, 2, 3, 13, "Winner + Love"),
               getStdPostCI(winner_angry_margin, 2, 3, 13, "Winner + Angry"),
               getStdPostCI(loser_love_margin, 2, 3, 13, "Loser + Love"),
               getStdPostCI(loser_angry_margin, 2, 3, 13, "Loser + Angry"),
               getStdPostCI(winner_love_poll, 2, 3, 13, "Winner + Love"),
               getStdPostCI(winner_angry_poll, 2, 3, 13, "Winner + Angry"),
               getStdPostCI(loser_love_poll, 2, 3, 13, "Loser + Love"),
               getStdPostCI(loser_angry_poll, 2, 3, 13, "Loser + Angry"),
               getStdPostCI(winner_love_surprise, 2, 3, 13, "Winner + Love"),
               getStdPostCI(winner_angry_surprise, 2, 3, 13, "Winner + Angry"),
               getStdPostCI(loser_love_surprise, 2, 3, 13, "Loser + Love"),
               getStdPostCI(loser_angry_surprise, 2, 3, 13, "Loser + Angry"))

citab$intr <- rep(c("(a) Post Election × Margin of Victory", 
                    "(b) Post Election × Poll Prediction",
                    "(c) Post Election × Surprise"),
                  each=4)

citab$label <- factor(citab$label, levels=rev(unique(citab$label)))
citab$intr <- factor(citab$intr, levels=unique(citab$intr))

# figure A.1
ggplot(data = citab, 
       aes(x = est, y = label)) + 
  facet_wrap(. ~ intr, ncol = 1) +
  geom_errorbar(data = citab, aes(xmin = lwr, xmax = upr),
                lwd=0.8, width = 0.1) +
  geom_vline(xintercept=0, lwd=0.5, lty=2) +
  geom_point(cex=4, pch=19) +
  xlab("Posterior Mean and 95% Credible Interval") +
  ylab("") +
  theme_bw() +
  theme(axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        axis.text=element_text(size=12),
        panel.spacing = unit(1, "lines"))

#ggsave("close_surprise.pdf", width=8, height=10)


# plot marginal effects
getME <- function(model, variable, label){
  posteriors <- posterior_samples(model, "b_post")

  est <- data.frame(t(sapply(1:length(variable), function(x){
    post_est <- unlist(posteriors[1] + posteriors[2] * variable[x])
    return(c(variable[x], 
             mean(post_est),
             quantile(post_est, 0.025),
             quantile(post_est, 0.975)))
  })))
  
  colnames(est) <- c("x", "est", "lwr", "upr")
  
  est$label <- label
  
  return(est)
}


# margin of victory
round(range(president_df$marginvictory_std, na.rm=TRUE))
mv_sim <- seq(-2, 41, length=100)

metab1 <- rbind(getME(winner_love_margin, mv_sim, "(a) Winner + Love"),
                getME(winner_angry_margin, mv_sim, "(b) Winner + Angry"),
                getME(loser_love_margin, mv_sim, "(c) Loser + Love"),
                getME(loser_angry_margin, mv_sim, "(d) Loser + Angry"))

metab1$label <- factor(metab1$label, levels=unique(metab1$label))

# figure A.2
ggplot(data = metab1, 
       aes(x = x, y = est)) + 
  geom_ribbon(data = metab1, aes(ymin = lwr, ymax = upr), alpha=0.3) +
  geom_hline(yintercept=0, lwd=0.5, lty=2) +
  geom_line(lwd=1.5) +
  facet_wrap(. ~ label, ncol = 4) +
  xlab("Margin of Victory (Standardized)") +
  ylab("Marginal Effect of Post Election") +
  theme_bw() +
  theme(axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15))

#ggsave("marginaleffect_margin.pdf", width=12, height=6)


# opinion poll prediction
round(range(president_df$poll_predict_std, na.rm=TRUE))
poll_sim <- seq(-5, 42, length=100)

metab2 <- rbind(getME(winner_love_poll, poll_sim, "(a) Winner + Love"),
                getME(winner_angry_poll, poll_sim, "(b) Winner + Angry"),
                getME(loser_love_poll, poll_sim, "(c) Loser + Love"),
                getME(loser_angry_poll, poll_sim, "(d) Loser + Angry"))

metab2$label <- factor(metab2$label, levels=unique(metab2$label))

# figure A.3
ggplot(data = metab2, 
       aes(x = x, y = est)) + 
  geom_ribbon(data = metab2, aes(ymin = lwr, ymax = upr), alpha=0.3) +
  geom_hline(yintercept=0, lwd=0.5, lty=2) +
  geom_line(lwd=1.5) +
  facet_wrap(. ~ label, ncol = 4) +
  xlab("Poll Prediction (Standardized)") +
  ylab("Marginal Effect of Post Election") +
  theme_bw() +
  theme(axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15))

#ggsave("marginaleffect_poll.pdf", width=12, height=6)


# surprise
round(range(president_df$surprise, na.rm=TRUE))
surprise_sim <- seq(-22, 26, length=100)

metab3 <- rbind(getME(winner_love_surprise, surprise_sim, "(a) Winner + Love"),
                getME(winner_angry_surprise, surprise_sim, "(b) Winner + Angry"),
                getME(loser_love_surprise, surprise_sim, "(c) Loser + Love"),
                getME(loser_angry_surprise, surprise_sim, "(d) Loser + Angry"))

metab3$label <- factor(metab3$label, levels=unique(metab3$label))

# figure A.4
ggplot(data = metab3, 
       aes(x = x, y = est)) + 
  geom_ribbon(data = metab3, aes(ymin = lwr, ymax = upr), alpha=0.3) +
  geom_hline(yintercept=0, lwd=0.5, lty=2) +
  geom_line(lwd=1.5) +
  facet_wrap(. ~ label, ncol = 4) +
  xlab("Surprise") +
  ylab("Marginal Effect of Election") +
  theme_bw() +
  theme(axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15))

#ggsave("marginaleffect_surprise.pdf", width=12, height=6)
